#This to obtain:
#Figure_A8_a.jpg
#Figure_A8_b.jpg
#Figure_A8_c.jpg
#Figure_A8_d.jpg
#Figure_A8_e.jpg
#Figure_A8_f.jpg

#table_A8.tex
#table_A9.tex

library("broom")
library("dplyr")
library("grid")
library("tidyr")
library("sandwich")
library("lfe")
library("broom")
library("dplyr")
library("tidyr")
library("grid")
library("stargazer")
library("lmtest")
library("multiwayvcov")
library("ggplot2")
library("glue")
library("spdep")
library("stringr")
library("lemon")
library("readxl")
library("spatialreg")

#setwd("C:/Users/bogdanp/")
setwd("/Users/bgpopescu/")

#Reading Polygons
x<-st_layers(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb")


HRV_adm0 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm0_wgs")

HRV_adm3 <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                    layer="HRV_adm3")
HRV_adm3<-st_simplify(HRV_adm3, dTolerance = 0.005)
HRV_adm3<-subset(HRV_adm3, select = c(ID_2))
HRV_adm3_data = read_excel("./Dropbox/Legacies_Project/Analysis/data/merge.xlsx", sheet=1, col_names = TRUE, skip = 0)


final_sp<-left_join(HRV_adm3, HRV_adm3_data, by = c("ID_2"="ID_2"))
final_sp<-subset(final_sp, bosnia_border_mun==0)

#################################
#Function for spatial regression#
#################################
spatial_regression <- function(data_sample, y, x, unique_observation_id) {
  clean_x = clean_list(x, "quad")
  clean_x = clean_list(clean_x, "bfe")
  
  specification = as.formula(paste(glue("{y}~"), 
                                   paste(clean_x, collapse = "+")))
  bw<-subset(data_sample, select =c(x, y, unique_observation_id))
  #Removing NAs
  bw<-bw[!sf::st_is_empty(bw), ]
  bw<-na.omit(bw)
  #Step1: the function "poly2nb" builds a neighbours list based on
  list_nb=poly2nb(bw, row.names = dplyr::pull(bw, unique_observation_id))
  #Step2: Convert nb to listw type
  list_nb_w=tryCatch(nb2listw(list_nb, zero.policy=TRUE), error=function(err) NA)
  #Step3: Spatial Filtering
  reg1_res_sp<-tryCatch(SpatialFiltering(specification,
                                         data = bw,
                                         nb = list_nb,
                                         style="W",
                                         ExactEV = F,
                                         zero.policy = T,
                                         na.action=na.exclude),
                        error=function(err) NA)

  print("Done Spatial Filtering")
  #Step4: Include fitted results"
  bw_fitted<-tryCatch(cbind(bw, as.data.frame(fitted(reg1_res_sp)),
                            id = row.names(bw)), error=function(err) 
                            {print("NO NEIGHBORS"); data_sample})
  #Keeping track of indices
  #This will allow us to see exactly which observations were used in the regressions
  row.names(bw_fitted) <- bw_fitted$id
  bw_fitted<-subset(bw_fitted, select = -c(id))
  list_fitted_ei<-names(bw_fitted[, grepl("vec", names(bw_fitted))])
  list_fitted_ei<-list_fitted_ei[list_fitted_ei != "Shape"]
  all_vars<-append(x, list_fitted_ei)
  specification_with_eigenv = as.formula(paste(glue("{y}~"), 
                                               paste(all_vars, collapse = "+")))
  
  print(specification_with_eigenv)
  model1 <- tryCatch(lm(specification_with_eigenv, data = bw_fitted, na.action = na.exclude), error=function(err) NA)
  
  reg_moran<-moran.test(x=resid(model1),listw=list_nb_w, zero.policy = T)$statistic
  reg_moran_p<-moran.test(x=resid(model1),listw=list_nb_w, zero.policy = T)$p.value
  reg_moran_star<-paste(round(reg_moran, digits = 3), "",
                        ifelse(reg_moran_p<.01,"***",
                               ifelse(reg_moran_p<.05,"**",
                                      ifelse(reg_moran_p<.1,"*",
                                             ""))), "", "", sep = "")
  
  
  return(list(model1, reg_moran_star))
}


######################################
#Function for cleaning the list of FE#
#####################################

clean_list <- function(old_list, particle) {
  term <- paste(glue("^{particle}"))
  new_list = list() 
  for (i in old_list) {
    #print(i)
    temp_list = c()
    
    for (j in i){
      if (!grepl(term, j)){
        temp_list <- c(temp_list, j)
      }
    }
    l = length(new_list)
    new_list[[l+1]] <- temp_list
  }
  return(new_list)}

#clean_list(specifications, "SN_")


###############################
#Function for rerunning models#
###############################
rerun_regression<-function(regression_name, string_dv, orig_data_for_clusters, cluster_name){
  #Process:
  #-takes already existing regression objects, associates them with an 
  #"official" label, calculates robust clustered SE and CI based on a cluster
  #that is in the original dataset
  #and puts the results in a tibble
  #-this will be used for ggplot later
  
  #Input:
  #-Regression model object after lm
  #-The label for the DV that you want displayed in ggplot
  #-original dataset that contains thhe names of the clusters
  
  #Output:
  #-A tibble with regression coefficients from the model that you ran with:
  #--robust SE
  #--robust CI
  
  #Recreate dataframe from the regression object
  remake_data<-regression_name$model
  #Scale the DV
  remake_data[1]<-scale(remake_data[1])
  #Take the string of the DV and IV
  dv<-regression_name$terms[[2]]
  ivs<-regression_name$terms[[3]]
  ivs_str<-paste(ivs)[2]
  dv_str<-paste(dv, "~", collapse = "")
  dv_ivs_str<-paste(dv_str, ivs_str, sep=" ")
  #Create specification and run model
  specification = as.formula(dv_ivs_str)
  model <- lm(specification, data = remake_data)
  no_obs<-nobs(model)
  no_obs<-paste("Obs: ", as.character(no_obs), sep = "")
  #Calculate robust SE
  model_se<-se_robust_cluster(model, orig_data_for_clusters, cluster_name)
  #Create tibble
  model_tidy<-tidy(model)
  model_tidy$no_obs<-no_obs
  model_se_tidy<-tidy(model_se)
  names(model_se_tidy)[1] <- "term"
  names(model_se_tidy)[2] <- "robust.se"
  model_tidy2<-left_join(model_tidy, model_se_tidy, by = c("term"="term"))
  
  #Calculate CI
  df_ci<-as.data.frame(ci_robust_cluster(model, orig_data_for_clusters, cluster_name))
  df_ci <- cbind(term = rownames(df_ci), df_ci)
  names(df_ci)[2] <- "lb_95"
  names(df_ci)[3] <- "ub_95"
  
  #Merge CI to the dataframe
  model_tidy3<-left_join(model_tidy2, df_ci, by = c("term"="term"))
  #Only keep the treat variable
  model_tidy4<-subset(model_tidy3, term=="treat")
  model_tidy4$dv_official<-string_dv
  return(model_tidy4)
}

se_robust_cluster <- function(x, data_sample_shape, name_cluster){
  data_model<-x$model
  data_model$in_regression<-1
  data_model<-subset(data_model, select=c(in_regression))
  #Create index
  data_model <- tibble::rowid_to_column(data_model, "index")
  data2<-tibble::rowid_to_column(data_sample_shape, "index")
  #Merge by rowname or index
  data_merged <- left_join(data_model, data2, by=c("index" = "index"))
  #Remove any missing
  data_subset<-subset(data_merged, !is.na(in_regression))
  dat_clusters<-data_subset[[name_cluster]]
  res<-coeftest(x, vcov = cluster.vcov(x, dat_clusters, force_posdef=T))[, 2]
  return(res)}

ci_robust_cluster <- function(x, data_sample_shape, name_cluster){
  data_model<-x$model
  data_model$in_regression<-1
  data_model<-subset(data_model, select=c(in_regression))
  #Create index
  data_model <- tibble::rowid_to_column(data_model, "index")
  data2<-tibble::rowid_to_column(data_sample_shape, "index")
  #Merge by rowname or index
  data_merged <- left_join(data_model, data2, by=c("index" = "index"))
  #Remove any missing
  data_subset<-subset(data_merged, !is.na(in_regression))
  dat_clusters<-data_subset[[name_cluster]]
  res_ci<-coefci(x, vcov = cluster.vcov(x, dat_clusters, force_posdef=T))
  return(res_ci)}


#Analyzing data
data<-final_sp

#Step1: Defining main IVs
specifications <- list(c('treat', 
                         'POINT_X', 'POINT_Y', 
                         'zagreb_distance', 
                         'bfe1', 'bfe2', 'bfe3', 'bfe4', 'bfe5', 'bfe6', 'bfe7',
                         'quad1', 'quad2', 'quad3'),
                       c('treat', 
                         'river_density', 'tmp_mean_2010', 'trade_route_1450_density', 
                         'POINT_X', 'POINT_Y', 
                         'zagreb_distance', 
                         'bfe1', 'bfe2', 'bfe3', 'bfe4', 'bfe5', 'bfe6', 'bfe7',
                         'quad1', 'quad2', 'quad3'),
                       c('treat', 'krajna6_distance',
                         'zagreb_distance', 
                         'bfe1', 'bfe2', 'bfe3', 'bfe4', 'bfe5', 'bfe6', 'bfe7',
                         'quad1', 'quad2', 'quad3'),
                       c('treat', 'krajna6_distance',
                         'POINT_X', 'POINT_Y', 
                         'zagreb_distance', 
                         'bfe1', 'bfe2', 'bfe3', 'bfe4', 'bfe5', 'bfe6', 'bfe7',
                         'quad1', 'quad2', 'quad3'))

####################
#Running the Models#
####################

latlon<-spatial_regression(data, "pct_no_water", specifications[[1]], "ID_2")[[1]]
latlon_m<-spatial_regression(data, "pct_no_water", specifications[[1]], "ID_2")[[2]]

latlon_unbal<-spatial_regression(data, "pct_no_water", specifications[[2]], "ID_2")[[1]]
latlon_unbal_m<-spatial_regression(data, "pct_no_water", specifications[[2]], "ID_2")[[2]]

dist<-spatial_regression(data, "pct_no_water", specifications[[3]], "ID_2")[[1]]
dist_m<-spatial_regression(data, "pct_no_water", specifications[[3]], "ID_2")[[2]]

latlon_dist<-spatial_regression(data, "pct_no_water", specifications[[4]], "ID_2")[[1]]
latlon_dist_m<-spatial_regression(data, "pct_no_water", specifications[[4]], "ID_2")[[2]]

sample_latlon_dist_exc<-subset(data,zagreb_area==0)
latlon_dist_exc<-spatial_regression(sample_latlon_dist_exc, "pct_no_water", specifications[[1]], "ID_2")[[1]]
latlon_dist_exc_m<-spatial_regression(sample_latlon_dist_exc, "pct_no_water", specifications[[1]], "ID_2")[[2]]

sample_ex_riv<-subset(data, krajna6_river_NEAR_FID==6)
ex_riv<-spatial_regression(sample_ex_riv, "pct_no_water", specifications[[1]], "ID_2")[[1]]
ex_riv_m<-spatial_regression(sample_ex_riv, "pct_no_water", specifications[[1]], "ID_2")[[2]]

marketplaces<-spatial_regression(data, "marketplaces", specifications[[1]], "ID_2")[[1]]
marketplaces_m<-spatial_regression(data, "marketplaces", specifications[[1]], "ID_2")[[2]]

rat_lawyers_notaries_1857<-spatial_regression(data, "rat_lawyers_notaries_1857", specifications[[1]], "ID_2")[[1]]
rat_lawyers_notaries_1857_m<-spatial_regression(data, "rat_lawyers_notaries_1857", specifications[[1]], "ID_2")[[2]]

rat_medical_personel_1857<-spatial_regression(data, "rat_medical_personel_1857", specifications[[1]], "ID_2")[[1]]
rat_medical_personel_1857_m<-spatial_regression(data, "rat_medical_personel_1857", specifications[[1]], "ID_2")[[2]]

#Creating Samples for Maps
sample_latlon_dist_exc_yes<-subset(sample_latlon_dist_exc, treat==1)
sample_latlon_dist_exc_no<-subset(sample_latlon_dist_exc, treat==0)
sample_latlon_dist_exc_yes = st_transform(sample_latlon_dist_exc_yes, st_crs(HRV_adm0))
sample_latlon_dist_exc_no = st_transform(sample_latlon_dist_exc_no, st_crs(HRV_adm0))
mil_col_short_polygon <- st_read(dsn="./Dropbox/Legacies_Project/Analysis/data/data.gdb",
                                 layer="krajna_line_GCS")
mil_col_short_polygon = st_transform(mil_col_short_polygon, st_crs(HRV_adm0))
final_sp_incl_brd<-st_simplify(data, dTolerance = 500)
final_sp_incl_brd = st_transform(data, st_crs(HRV_adm0))


cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Civilian Territory" = alpha("blue", 0.5),
              "Excluded" = alpha("gray", 0.01))

map<-ggplot() + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = sample_latlon_dist_exc_yes, aes(fill="Military Territory"),
          linewidth = 0.1)+
  geom_sf(data = sample_latlon_dist_exc_no, aes(fill="Civilian Territory"),
          linewidth = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  ggtitle("(5) Excluding Zagreb")+
  labs(x = "Longitude", y="Latitude")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
    theme(axis.title=element_text(size=14),
          plot.title = element_text(hjust = 0.5),
          legend.position = c(1, 0),
          #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
          #and c(1,1) corresponds to the "top right" position.
          legend.box.background = element_rect(fill='white'),
          legend.background = element_blank(),
          legend.text=element_text(size=12),
          plot.margin = unit(c(0, 0, 0, 0), "null"),
          panel.spacing = unit(c(0, 0, 0, 0), "cm"))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/Figure_A8_a.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)


#Creating samples that excludes districts whose border coincides with borders
sample_excl_river<-subset(data, krajna6_river_NEAR_FID==6)
sample_excl_river_yes<-subset(sample_excl_river, treat==1)
sample_excl_river_no<-subset(sample_excl_river, treat==0)
sample_excl_river_yes = st_transform(sample_excl_river_yes, st_crs(HRV_adm0))
sample_excl_river_no = st_transform(sample_excl_river_no, st_crs(HRV_adm0))
mil_col_short_polygon = st_transform(mil_col_short_polygon, st_crs(HRV_adm0))


cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Civilian Territory" = alpha("blue", 0.5),
              "Excluded" = alpha("gray", 0.01))


map<-ggplot() + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = sample_excl_river_yes, aes(fill="Military Territory"),
          size = 0.1)+
  geom_sf(data = sample_excl_river_no, aes(fill="Civilian Territory"),
          size = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  ggtitle("(6) Excluding River-Borders")+
  labs(x = "Longitude", y="Latitude")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/Figure_A8_b.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)


#Remove border area which neighbors Bosnia
sample_non_border_area <- subset(data, bosnia_border_area == 0)
non_border_area<-spatial_regression(sample_non_border_area, "pct_no_water", specifications[[1]], "ID_2")[[1]]
non_border_area_m<-spatial_regression(sample_non_border_area, "pct_no_water", specifications[[1]], "ID_2")[[2]]

#Creating Samples for Maps
sample_non_border_area_yes<-subset(sample_non_border_area, treat==1)
sample_non_border_area_no<-subset(sample_non_border_area, treat==0)
sample_non_border_area_yes = st_transform(sample_non_border_area_yes, st_crs(HRV_adm0))
sample_non_border_area_no = st_transform(sample_non_border_area_no, st_crs(HRV_adm0))
mil_col_short_polygon = st_transform(mil_col_short_polygon, st_crs(HRV_adm0))


cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Civilian Territory" = alpha("blue", 0.5),
              "Excluded" = alpha("gray", 0.01))

map<-ggplot() + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = sample_non_border_area_yes, aes(fill="Military Territory"),
          size = 0.1)+
  geom_sf(data = sample_non_border_area_no, aes(fill="Civilian Territory"),
          size = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  ggtitle("(7) Excluding Border Regions")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/Figure_A8_c.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)


#Selecting areas which were not attacked by the Ottomans
sample_excl_ottoman_attacks_res<-subset(data, ottoman_attack==0)
excl_ottoman_attacks_res<-spatial_regression(sample_excl_ottoman_attacks_res, "pct_no_water", specifications[[1]], "ID_2")[[1]]
excl_ottoman_attacks_res_m<-spatial_regression(sample_excl_ottoman_attacks_res, "pct_no_water", specifications[[1]], "ID_2")[[2]]
summary(excl_ottoman_attacks_res)

#Creating Samples for Maps
sample_excl_ottoman_attacks_res_yes<-subset(sample_excl_ottoman_attacks_res, treat==1)
sample_excl_ottoman_attacks_res_no<-subset(sample_excl_ottoman_attacks_res, treat==0)
sample_excl_ottoman_attacks_res_yes = st_transform(sample_excl_ottoman_attacks_res_yes, st_crs(HRV_adm0))
sample_excl_ottoman_attacks_res_no = st_transform(sample_excl_ottoman_attacks_res_no, st_crs(HRV_adm0))
mil_col_short_polygon = st_transform(mil_col_short_polygon, st_crs(HRV_adm0))


cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Civilian Territory" = alpha("blue", 0.5),
              "Excluded" = alpha("gray", 0.01))

map<-ggplot() + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = sample_excl_ottoman_attacks_res_yes, aes(fill="Military Territory"),
          linewidth = 0.1)+
  geom_sf(data = sample_excl_ottoman_attacks_res_no, aes(fill="Civilian Territory"),
          linewidth = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  ggtitle("(8) Excluding Ottoman attacked Regions")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/Figure_A8_d.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)

#Removing mushroom area
sample_border_area<-subset(data, mushroom==0)
border_area<-spatial_regression(sample_border_area, "pct_no_water", specifications[[1]], "ID_2")[[1]]
border_area_m<-spatial_regression(sample_border_area, "pct_no_water", specifications[[1]], "ID_2")[[2]]

#Creating Samples for Maps
sample_border_area_yes<-subset(sample_border_area, treat==1)
sample_border_area_no<-subset(sample_border_area, treat==0)
sample_border_area_yes = st_transform(sample_border_area_yes, st_crs(HRV_adm0))
sample_border_area_no = st_transform(sample_border_area_no, st_crs(HRV_adm0))
mil_col_short_polygon = st_transform(mil_col_short_polygon, st_crs(HRV_adm0))


cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Civilian Territory" = alpha("blue", 0.5),
              "Excluded" = alpha("gray", 0.01))

map<-ggplot() + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = sample_border_area_yes, aes(fill="Military Territory"),
          linewidth = 0.1)+
  geom_sf(data = sample_border_area_no, aes(fill="Civilian Territory"),
          linewidth = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  ggtitle("(9) Only Border Regions")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/Figure_A8_e.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)

#Removing area which was not populated by the Orthodox
sample_immigration_type_orthodox<-subset(data,immigration_type_orthodox==0)
immigration_type_orthodox<-spatial_regression(sample_immigration_type_orthodox, "pct_no_water", specifications[[1]], "ID_2")[[1]]
immigration_type_orthodox_m<-spatial_regression(sample_immigration_type_orthodox, "pct_no_water", specifications[[1]], "ID_2")[[2]]

#Creating Samples for Maps
sample_immigration_type_orthodox_yes<-subset(sample_immigration_type_orthodox, treat==1)
sample_immigration_type_orthodox_no<-subset(sample_immigration_type_orthodox, treat==0)
sample_immigration_type_orthodox_yes = st_transform(sample_immigration_type_orthodox_yes, st_crs(HRV_adm0))
sample_immigration_type_orthodox_no = st_transform(sample_immigration_type_orthodox_no, st_crs(HRV_adm0))
mil_col_short_polygon = st_transform(mil_col_short_polygon, st_crs(HRV_adm0))


cols.fill = c("Military Territory" = alpha("red", 0.5), 
              "Civilian Territory" = alpha("blue", 0.5),
              "Excluded" = alpha("gray", 0.01))

map<-ggplot() + 
  geom_sf(data = HRV_adm0, fill = NA, colour = "grey10", linewidth = 0.06)+
  geom_sf(data = final_sp_incl_brd, aes(fill="Excluded"), colour = "grey60", linewidth = 0.06)+
  geom_sf(data = sample_immigration_type_orthodox_yes, aes(fill="Military Territory"),
          linewidth = 0.1)+
  geom_sf(data = sample_immigration_type_orthodox_no, aes(fill="Civilian Territory"),
          linewidth = 0.1)+
  geom_sf(data = mil_col_short_polygon, colour = "blue", fill = NA, linewidth = 0.5)+
  scale_fill_manual(name = "Sample", values = cols.fill)+
  theme_bw()+
  labs(x = "Longitude", y="Latitude")+
  ggtitle("(10) Excluding Orthodox Regions")+
  coord_sf(ylim=c(42.4,46.669853), xlim=c(13.556495, 19.447201))+
  theme(plot.title = element_text(hjust = 0.5))+
  #scale_fill_manual("", values = my.cols, guide = "legend") +
  theme(axis.title=element_text(size=14),
    plot.title = element_text(hjust = 0.5),
    legend.position = c(1, 0),
    #Legend.position values should be between 0 and 1. c(0,0) corresponds to the "bottom left"
    #and c(1,1) corresponds to the "top right" position.
    legend.box.background = element_rect(fill='white'),
    legend.background = element_blank(),
    legend.text=element_text(size=12))
map2<-map+ggsn::scalebar(x.min = 18.5, x.max = 19.1,
                         #Above you mention how long the scalebar should be
                         y.min = 46.4, y.max = 46.7,
                         #Above you mention how thick the scalebar should be
                         dist = 50, dist_unit = "km",
                         transform = T, model = "WGS84",
                         location = "topright",
                         st.size = 3,
                         #Above you have the font size of the numbers below the scalebar
                         st.dist =0.4,
                         #Above you have the distance between the bar and the text, as a proportion of the y axis.
                         height=0.2)
# #Above you have anumber between 0 and 1 to indicate the height 
# #of the scale bar, as a proportion of the y axis
map3<-reposition_legend(map2, 'bottom left')

ggsave(map3,file="./Dropbox/Legacies_Project/Paper/figures/samples/Figure_A8_f.jpg", 
       height=15, width=15, 
       units = "cm", dpi=300)


mean_pct_no_water<-round(mean(data$pct_no_water, na.rm=T), digits = 3)
sd_pct_no_water<-round(sd(data$pct_no_water, na.rm=T), digits = 3)


column_labels= c("\\shortstack{Lat-Lon}", 
                 "\\shortstack{Lat-Lon +\\\\ Controls}",
                 "\\shortstack{Distance\\\\ Brd.}",
                 "\\shortstack{Lat-Lon and\\\\ Distance Bnd.}",
                 "\\shortstack{Excluding\\\\ Zagreb}",
                 "\\shortstack{Excluding\\\\ River-Borders}",
                 "\\shortstack{Excluding\\\\ Border\\\\ Regions}",
                 "\\shortstack{Excluding\\\\ Ottoman attacked\\\\ Regions}",
                 "\\shortstack{Only\\\\ Border\\\\ Regions}",
                 "\\shortstack{Excluding \\\\ Orthodox\\\\ Regions}")

list_models<-list(latlon,
               latlon_unbal,
               dist,
               latlon_dist,
               latlon_dist_exc,
               ex_riv,
               non_border_area,
               excl_ottoman_attacks_res,
               border_area,
               immigration_type_orthodox)

pct_no_water <-stargazer(list_models,
                        column.labels= column_labels,
                        keep = c("treat", 'river_density', 'tmp_mean_2010', 'trade_route_1450_density'),
                        se = list(se_robust_cluster(latlon, data, "NAME_2"),
                                  se_robust_cluster(latlon_unbal, data, "NAME_2"),
                                  se_robust_cluster(dist, data, "NAME_2"),
                                  se_robust_cluster(latlon_dist, data, "NAME_2"),
                                  se_robust_cluster(latlon_dist_exc, sample_latlon_dist_exc, "NAME_2"),
                                  se_robust_cluster(ex_riv, sample_excl_river, "NAME_2"),
                                  se_robust_cluster(non_border_area, sample_non_border_area, "NAME_2"),
                                  se_robust_cluster(excl_ottoman_attacks_res, sample_excl_ottoman_attacks_res, "NAME_2"),
                                  se_robust_cluster(border_area, sample_border_area, "NAME_2"),
                                  se_robust_cluster(immigration_type_orthodox, sample_immigration_type_orthodox, "NAME_2")),
                        covariate.labels=c("Military Territory", "River Density", 
                                           "Temperature Avg.", "Trade Route Dens., 1450"),
                        dep.var.caption = "Dependent variable: Pct. Dwellings No Water, 2011.
                                Specification:",
                        dep.var.labels.include = F,
                        omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                        add.lines=list(c("Mean", rep(mean_pct_no_water, 10)),
                                       c("SD", rep(sd_pct_no_water, 10)),
                                       c("Boundary FE", rep("Yes", 10)),
                                       c("Region FE", rep("Yes", 10)),
                                       c("Moran eigenvectors", rep("Yes", 10)),
                                       c("Moran's I Residual", latlon_m,
                                         latlon_unbal_m,
                                         dist_m,
                                         latlon_dist_m,
                                         latlon_dist_exc_m,
                                         ex_riv_m,
                                         non_border_area_m,
                                         excl_ottoman_attacks_res_m,
                                         border_area_m,
                                         immigration_type_orthodox_m)))
note_latex <- "\\multicolumn{11}{l} {\\parbox[t]{30cm}{ \\footnotesize \\textit{Notes:}
Coefficients and county-clustered robust standard errors in parantheses from OLS regression.
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01. The unit of analysis is municipality (opcina).
Data on percentage access to water is from the Croatian Bureau of Statistics (2017, pp. 157-211),
based on data from the 2011 census. See codebook in the appendix for data sources and how the variables were calculated.}}"
pct_no_water [grepl("Note", pct_no_water)] <- note_latex
#cat(pct_no_water, file="pct_no_water.tex", sep="\n")
cat(pct_no_water, file="./Dropbox/Legacies_Project/Paper/tables/table_A8.tex", sep="\n")



mean_marketplaces<-round(mean(data$marketplaces, na.rm=T), digits = 3)
sd_marketplaces<-round(sd(data$marketplaces, na.rm=T), digits = 3)
mean_rat_lawyers_notaries_1857<-round(mean(data$rat_lawyers_notaries_1857, na.rm=T), digits = 3)
sd_rat_lawyers_notaries_1857<-round(sd(data$rat_lawyers_notaries_1857, na.rm=T), digits = 3)
mean_rat_medical_personel_1857<-round(mean(data$rat_medical_personel_1857, na.rm=T), digits = 3)
sd_rat_medical_personel_1857<-round(sd(data$rat_medical_personel_1857, na.rm=T), digits = 3)


column_labels= c("\\shortstack{No. Marketplaces, 1857}", 
                 "\\shortstack{Prop. Lawyers and Notaries, 1857}",
                 "\\shortstack{Prop. Doctors, 1857}")

table_A9 =stargazer(marketplaces,
                   rat_lawyers_notaries_1857,
                   rat_medical_personel_1857,
                   column.labels= column_labels,
                   keep = c("treat"),
                   se = list(se_robust_cluster(marketplaces, data, "NAME_2"),
                             se_robust_cluster(rat_lawyers_notaries_1857, data, "NAME_2"),
                             se_robust_cluster(rat_medical_personel_1857, data, "NAME_2")),
                   covariate.labels=c("Military Territory"),
                   dep.var.caption = "Dependent variable:",
                   dep.var.labels.include = F,
                   omit.stat=c("LL","ser","f", "rsq"), no.space=TRUE, float=FALSE,
                   add.lines=list(c("Mean", mean_marketplaces, mean_rat_lawyers_notaries_1857, mean_rat_medical_personel_1857),
                                  c("SD", 
                                    sd_marketplaces,
                                    sd_rat_lawyers_notaries_1857,
                                    sd_rat_medical_personel_1857),
                                  c("Boundary FE", rep("Yes", 10)),
                                  c("Region FE", rep("Yes", 10)),
                                  c("Moran eigenvectors", rep("Yes", 10)),
                                  c("Moran's I Residual", marketplaces_m,
                                    rat_lawyers_notaries_1857_m,
                                    rat_medical_personel_1857_m)))
note_latex <- "\\multicolumn{4}{l} {\\parbox[t]{18cm}{ \\footnotesize \\textit{Notes:}
Coefficients and county-clustered robust standard errors in parantheses from OLS regression.
$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01. The unit of analysis is municipality (opcina) after overlaying
them with 1857 units (bezirke). The data for the dependent variables is from Austrian Ministry of Internal Affairs (1859). 
Dependent variables are standardized.
See codebook for details.}}"
table_A9 [grepl("Note", table_A9)] <- note_latex
cat(table_A9, file="./Dropbox/Legacies_Project/Paper/tables/table_A9.tex", sep="\n")
